TIME  ACCURATE  COMPUTATION  OF  UNSTEADY  SHOCK  TUNNEL 
FLOW  WITH  COUPLED  DIAPHRAGM  RUPTURE  MECHANICS 


D.  Scott  McRae 
M.  A.  Zikry 

Department  of  Mechanical  and  Aerospace  Engineering 
North  Carolina  State  University 
Raleigh,  North  Carolina  27695-7910 


29  October  1 999 

Final  Technical  Report  for  the  Period  1  August  1995  -  31  July  1999 
AFOSR  Grant  F49620-95- 1-0374 


Distribution  Unlimited 


Prepared  for 


Air  Force  Office  of  Scientific  Research 
AFOSR/MN 

801  N.  Randolph  Street,  Room  732 
Arlington,  VA  22203 


Dll 


20000608  101 


REPORT  DOCUMENTATION  PAGE 


Publk  reporting  burden  fw  tins  collection  ofinfomiation  is  estimated  to  average  1  hour  per  response,  including  the  tune  for  reviewing  in 

and  completing  and  reviewing  this  collection  of  infonnation.  Send  comments  regarding  this  burden  estimate  or  any  other  a^)ect  of  this  <  ^ 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1 204,  Arlington,  VA 
(0704-0188),  Washington,  DC  20503 _ _ _ - _ _ 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  ainu  - - - - 

29  OCTOBER  1999  FINAL  TECHNICAL  REPORT-  1  AUGUST  1995  TO  31  JULY  1999 


AFRL-SR-BL-TR-OO- 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS 

Time  Accurate  Computation  of  Unsteady  Shock  Tunnel  Flow  with  Coupled  AFOSR  Grant  F49620-95-I-0374 

Diaphragm  Rupture  Mechanics 


needed, 
ashington 
tion  Project 


2.  REPORT  DATE 
29  OCTOBER  1999 


6.  AUTHOR(S) 

D.  S.  McRae 
M.  A.  Zikry 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

N.C.  state  University 
Department  of  Mechanical  and 
Aerospace  Engineering. 

Box  7910 

Raleiah,  NC  27695-791 0  _ _ 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Dr.  Len  Sakell 
AFOSR/MN 

Air  Force  Office  of  Scientific  Research 
801  N.  Randolph  Street,  Room  732 

Ariinoton.  VA  22203  _ 


11.  SUPPLE.MENTARY  NOTES 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBITION  /  AVAILABILITY  STATEMENT 

Unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTR/VCT  (Maximtun  200  Words) 

This  report  reviews  work  performed  toward  the  goal  of  simulating  time  accurately  the  initiation  of  flow  in  a 
hypersonic  shock  tunnel,  including  coupled  diaphragm  rupture  mechanics.  A  solution-adaptive  mesh  Naver-Stokes  code 
was  modified  to  improve  mesh  distribution  for  rapidly  moving  flow  features,  to  include  boundary  surface  motion,  and  to 
incorporate  weight  function  and  stability  enhancements.  A  simulation  of  an  unadapted  grid  blade-type  shock  tube 
diaphragm  opening  actuated  by  algorithm  was  completed.  A  repeat  of  these  simulations  with  mesh  adaptation  revealed 
problems  with  the  adaptive  algorithm  related  to  the  interaction  with  the  grid  block  structure  that  were  only  partially 
alleviated  by  limiting  mesh  movement  at  the  blade  tip. 

Plane  strain  and  axisymmetric  finite-element  programs  have  been  developed  for  the  investigation  of  dynamic 
inelastic  deformation  and  rupture  of  a  shock-tube  diaphragm.  Results  indicate  that  crack  nucleation  occurs  when 
compressive  waves  at  high  pressures  reflect  off  the  rear  of  the  diaphragm  as  tensile  waves  that  are  a  precursor  to  crack 
initiation.  Results  from  this  study,  indicate  that  stress  fields,  stress  gradients,  and  plastic  strains  at  the  free-edge  or  at  the 
crack  or  notch  are  physically  large  enough,  at  the  higher  pressures,  to  result  in  crack  nucleation,  diaphragm  rupture,  and 
material  separation 


14.  SUBJECT  TERMS 

CFD,  shock  tunnel,  solid  mechanics,  inelastic  strain,  crack  propagation,  stress  waves 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 


17,  SECURITY  CLASSIFICATION 

18.  SECURITY 

OF  REPORT 

CLASSIFICATION 

Unclassified 

OF  THIS  PAGE  Unclassified 

15.  NUMBER  OF  PAGES 


lb.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 

None 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Presoribed  by  ANSI  Std.  Z39-18 
298-102 


Executive  Summary 


Technical  Review 

This  report  reviews  work  performed  toward  the  goal  of  simulating  time  accurately 
the  initiation  of  flow  in  a  hypersonic  shock  tunnel,  including  coupled  diaphragm  rupture 
mechanics.  Two  codes  are  under  development  to  investigate  coupled  flow  and 
diaphragm  simulation.  The  fluids  code  was  based  on  an  existing  solution-adaptive 
mesh  Naver-Stokes  code  that  had  been  applied  to  compute  the  unsteady  flow  over 
stalled  compressor  blades  arranged  in  a  cascade.  Initial  testing  of  this  code  revealed 
that  the  mesh  movement  algorithm,  which  had  performed  well  when  applied  to 
compressor  cascades  and  to  a  self  excited  oscillatory  flow,  was  inadequate  for  the 
large,  high-  velocity  feature  movement  present  in  shock  tube-  type  flows.  Modifications 
were  made  to  effect  improved  mesh  distribution  for  rapidly  moving  flow  features,  to 
include  boundary  surface  motion  in  both  solver  and  adapter,  and  to  incorporate  weight 
function  and  stability  enhancements  from  the  NCSU  adaptive  mesh  code  SIERRA 
(Laflin).  The  mesh  movement  and  block  boundary  modifications  were  verified  by  solving 
the  Riemann  problem  for  conditions  with  shock  speeds  similar  to  those  to  be 
encountered  in  an  actual  shock  tunnel.  The  boundary  movement  changes  were  verified 
by  a  piston  movement  computation,  whereas  both  of  these  computations  served  to 
verify  the  weight  function  and  stability  enhancements.  Next  a  simulation  of  a  shock  tube 
diaphragm  opening  without  grid  adaptation  was  completed.  A  blade-  type  diaphragm 
was  chosen  to  minimize  the  difficulties  associated  with  hanging  blocks  during 
adaptation.  When  the  diaphragm  was  actuated  by  algorithm,  these  simulations  revealed 
two  different  flow  initiation  mechanisms,  depending  on  the  ratio  of  surface  movement  to 
fluid  velocity.  Next  these  simulations  were  repeated  with  mesh  adaptation.  Problems 
were  encountered  related  to  the  interaction  of  the  grid  block  structure  with  the  adaptive 
algorithm.  Limiting  mesh  movement  in  the  vicinity  of  the  blade  tip  only  partially  alleviated 
these  problems.  Analysis  showed  that  one  of  (or  a  combination  of)  two  further 
developments  would  be  necessary  to  complete  the  fluids  code  development.  The  first 
was  to  allow  non-continuous  grids  across  block  boundaries.  This  was  needed  to  allow 
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each  block  to  undergo  independent  grid  deformation  without  imposing  the  same 
deformation  on  neighboring  blocks.  The  second  requirement  was  dictated  by  the  grid 
deformation  upon  full  opening  of  the  diaphragm.  An  adaptive  unstructured  grid  code 
would  permit  this  physical  deformation  while  maintaining  grid  viability.  Both 
developments  are  being  carried  fonward  under  other  research  projects,  and  will  be 
available  for  further  exploration  of  this  problem. 

During  the  solid  mechanics  research,  plane  strain  and  axisymmetric  finite- 
element  programs  have  been  developed  for  the  investigation  of  dynamic  inelastic 
deformation  and  rupture  of  a  shock-tube  diaphragm.  Modifications  have  also  been 
applied  for  the  characterization  of  damage  progression  and  failure  evolution  pertaining 
to  diaphragm  rupture.  The  dynamic  constitutive  formulations  that  were  developed 
account  for  high  pressures,  strain-rate,  and  thermal  effects.  Specialized  numerical 
techniques  have  been  developed  that  allow  tracking  of  propagating  stress  waves  in 
structural  components.  This  is  significant,  because  our  results  indicate  that  crack 
nucleation  occurs  when  compressive  waves  at  high  pressures  reflect  off  the  rear  of  the 
diaphragm  as  tensile  waves  that  are  a  precursor  to  crack  initiation.  This  two  dimensional 
wave  reflection  was  used  as  a  ductile  failure  criterion  to  signify  crack  initiation  for  high 
pressure  loading  conditions.  Using  this  failure  criterion,  stresses  that  correspond  to  the 
fracture  strength  of  the  material  have  been  predicted.  The  generalized  inelastic  rate- 
dependent  phenomenological  plasticity  formulation  that  accounts  for  noncoaxial  effects 
such  as  thermal  softening  and  kinematic  and  isotropic  hardening  has  been  used,  with 
the  newly  developed  failure  criteria,  to  investigate  damage  progression  and  failure 
evolution  in  rupturing  diaphragms.  Unnotched,  notched,  and  cracked  geometries  have 
been  used  to  analyze  the  inelastic  strain  fields  that  could  be  associated  with  either  the 
gradual  opening  (petal  opening),  or  the  instantaneous  rupturing  of  the  diaphragm.  New 
methods  have  developed  to  track  the  failure  due  to  the  reflected  waves.  Results  from 
this  study  indicate  that  stress  fields,  stress  gradients,  and  plastic  strains  at  the  free-edge 
or  at  the  crack  or  notch  are  physically  large  enough,  at  the  higher  pressures,  to  result  in 
crack  nucleation,  diaphragm  rupture,  and  material  separation.  This  buildup  in  the 
opening  mode  stress  is  directly  related  to  the  reflected  tensile  waves.  Based  on  the 
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newly  developed  solid  mechanics  code  for  structural  integrity,  the  onset  of  failure  was 
predicted  and  used  to  model  petal  opening. 

Graduate  Student  Retention 

Research  progress  in  this  project  was  limited  during  the  final  two  years  by  the 
very  issue  that  the  AASERT  program  was  intended  to  improve;  the  reduced  number  of 
U.S.  citizens  entering  graduate  school  in  the  core  engineering  disciplines.  Of  the  two 
students  funded  initially  by  this  grant,  the  first,  after  completing  the  MS  degree  and 
achieving  solid  progress  toward  developing  a  fluids  code  for  this  research,  chose  to 
accept  a  lucrative  software  related  position  with  E-Systems  rather  than  continue  for  the 
PhD.  The  second  student  who  performed  the  solid  mechanics  research,  after 
completing  an  MS,  chose  to  explore  other  research  projects  even  though  continuing  for 
the  PhD.  In  normal  times,  two  U.S.  students  would  have  been  recruited  from  our  CFD 
and  solid  mechanics  sequences  and  the  work  would  have  continued  with  some 
slowdown  for  training  and  code  familiarization.  However  no  students  from  these  groups 
were  immediately  available.  For  the  fluids  research,  a  student  was  found,  in  the  process 
of  returning  to  school  after  a  break  in  a  computational  astrophysics  MS  program,  who 
was  interested  in  a  dual  AE/Physics  degree  (BSAE  from  UVA).  However,  considerable 
training  in  CFD  was  needed  and  the  student  subsequently  took  a  full  time  software 
related  job  off  campus  because  of  financial  obligations  that  the  RA  stipend  would  not 
support.  Although  some  work  on  the  fluids  code  continued  beyond  this  point,  little  true 
progress  was  made  toward  the  final  goal.  In  the  case  of  the  solid  mechanics  student, 
he  returned  to  the  project  after  a  one  year  hiatus  and  continued  work  on  the  crack 
propagation  and  deformation  study  reported  herein.  However,  no  coupling  between  the 
two  codes  was  accomplished. 

Student  Retention  Recommendation 

Allowance  for  greatly  increased  graduate  student  stipends  should  be  included  in 
future  AFOSR/AASSERT  grants.  The  initial  reduction  in  undergraduate  student 
population  in  the  1990’s  was  clearly,  for  Aerospace  Engineering,  due  to  the  negative 
perceptions  of  job  availability  upon  graduation.  However,  now  that  the  undergraduate 
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AE  enrollment  is  rebounding,  the  increased  availability  of  lucrative  positions  at  both  the 
BS  and  MS  levels  has  kept  AE  graduate  enrollment  depressed.  It  is  unrealistic  to  expect 
many  qualified  students  to  forego  $30K  to  $50K  in  salary  differential  yearly  to  enter  or 
remain  in  graduate  school. 

Project  Personnel 

D.  Scott  McRae,  PI 
M.  A.  Zikry,  Co-PI 

Funded  Graduate  Research  Assistants 

Jared  Baucom,  Ph.D.  student  (Degree  expected  12/2000) 

Darren  White,  completed  M.S.  degree,  August  1997. 

G.  Christopher  Jenkins,  M.S.  student  replacing  D.  White,  August  1997(Degree  expected 
6/2000) 

Publications 

M.  A.  Zikry  and  J.  Baucom  (1997),  High  Pressure  Shear  Strain  Localization  in 
Structural  Steel,  in  Applications  of  Engineering  Plasticity,  pp.  371-376,  ed:  A.  Khan, 
Elsevier  Press.  A  presentation  was  also  given  at  the  Sixth  International  Plasticity 
Conference  held  in  Juneau,  Alaska,  July,  1997. 

J.  N.  Baucom  and  M.  A.  Zikry  (1999)  “Perturbation  analysis  of  high 
strain-rate  shear  localization  in  B.C.C.  crystalline  materials”  Acta 
Mechanica  137,  109-129. 

Thesis 

White,  D.  A.,  “Computation  and  Analysis  of  Hypersonic  Shock  Tubes  Including 
Diaphragm  Rupture  Effects,”  M.  S.  Thesis,  N.  C.  State  University,  August  1997. 
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FLUID  MECHANICS  SIMULATION 


This  section  will  review  work  performed  toward  modification  and  development  of 
the  multiblock  2-D  Navier-Stokes  that  was  continued  from  a  prior  research  project.  The 
primary  tasks  have  included  improvement  of  the  adaptive  algorithm  to  resolve  very 
rapidly  moving  solution  features,  installation  of  moving  surface  boundary  conditions  in 
both  the  solver  and  the  adapter,  incorporation  of  stability  and  accuracy  enhancements 
developed  by  Laflin  (SIERRA),  and  an  initial  simulation  of  an  actual  shock  tube 
diaphragm  opening  with  flow  initiation.  The  description  of  these  tasks  was  taken 
primarily  from  White.  Each  will  be  discussed  individually  and  results  will  be  included 
where  appropriate. 

Rapidly  Moving  Features  (Grid  Convergence) 

The  mass  weighted  algorithm  of  Eiseman  is  used  in  a  center  of  mass  approach 
in  the  NCSU  adaptive  algorithm.  Expansion  of  the  mass  weighted  algorithm  in  a  Taylor 
series  yields  an  elliptic  system  of  partial  differential  equations.  The  fundamental 
mathematical  character  of  an  elliptic  system  must  be  kept  in  mind  when  using  this 
approach  for  dynamic  grid  adaption.  In  previous  work  by  Benson,  Ingram,  and  Neaves, 
1  to  5  explicit  update  steps  were  been  used  to  solve  for  the  new  grid  point  coordinates. 
In  the  majority  of  these  cases  the  new  solution  was  only  a  small  perturbation  of  that  on 
the  previous  grid,  thereby  requiring  that  the  grid  node  move  only  a  small  distance  per 
time  step.  This  offers  the  benefit  of  being  much  more  computationally  affordable  than 
actually  driving  the  elliptical  PDE  system  to  convergence  at  each  time  step.  For  steady 
state  solutions  or  solutions  with  well  resolved  time  histories  (i.e.,  small  time  steps),  this 
approach  has  been  shown  to  produce  excellent  results  (Benson,  Ingram,  Neaves,  and 
Laflin).  However,  poor  results  were  obtained  when  the  standard  number  of  mass- 
weighted  algorithm  passes  were  used  in  the  shock  tube  calculations,  especially  when 
the  shock  and  expansion  waves  were  highly  resolved.  When  the  waves  moved  into  the 
evacuated  area  in  one  time  step,  the  resolution  was  not  sufficient  to  maintain  consistent 
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solution  accuracy.  After  some  analysis,  it  was  discovered  that  iterating  the  elliptical 
problem  until  convergence  was  approached  resulted  in  a  much  more  even  distribution 
of  mesh  nodes  in  regions  of  constant  weight  function. 

The  elliptic  center  of  mass  algorithm  does  provide  the  necessary  propagation  of 
information  throughout  the  domain  provided  that  the  solution  is  allowed  to  converge 
sufficiently.  For  full  mathematical  ellipticity,  the  effect  of  moving  one  grid  point  in  the 
domain  must  influence  the  position  of  every  other  grid  point  within  the  domain. 

Obviously  a  few  explicit  point-Jacoby  iteration  steps  cannot  provide  the  necessary 
global  propagation  of  solution  information. 

The  needed  convergence  of  the  grid  equations  can  be  achieved  through  an 
explicit  point  Jacoby  algorithm  only  if  a  sufficient  number  of  iterations  are  performed. 
For  increased  efficiency,  a  successive  over-relaxation  (SOR)  scheme  was  implemented. 
For  highly  dynamic  cases  requiring  many  iterations  to  relax  the  mesh,  an  implicit 
alternating  direction  line  SOR  scheme  was  used  to  further  speed  up  the  numerical 
propagation  of  information  throughout  the  domain.  Generally  speaking,  the  purpose  of 
the  explicit  SOR  and  implicit  algorithms  was  to  increase  the  numerical  domain  of 
dependence  of  the  solution,  not  to  provide  increased  local  point  movement.  If  the  flow 
features  being  tracked  are  moving  rapidly  enough,  an  explicit  SOR  or  an  implicit 
adaption  algorithm  can  result  in  locally  high  grid  velocities.  Special  care  must  therefore 
be  taken  with  these  schemes  to  limit  the  range  of  point  movement  to  ensure  that  grid 
lines  do  not  cross.  In  order  to  prevent  crossover,  a  grid  point  must  not  be  allowed  to 
move  outside  its  dual  cell  at  each  iteration  of  the  grid  solver.  A  numerical  experiment 
demonstrating  the  significance  of  grid  convergence  is  presented  in  the  Results  Section. 

Since  one  of  the  goals  of  this  research  was  to  couple  the  diaphragm  opening  with 
the  fluid  dynamics  using  a  solution  adaptive  mesh,  means  were  developed  to  adjust  for 
boundary  movement  while  adapting  the  mesh  internally.  The  general  integral  governing 
equation  for  fluids  can  be  written 


where  jc  is  a  general  vector  denoting  the  movement  of  a  specific  cell  side  relative  to  the 
inertial  mesh.  In  the  original  adaptive  algorithm,  this  vector  was  non-zero  only  in 
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response  to  changes  in  the  solution.  However,  when  the  outer  boundary  of  the  mesh 
must  be  moved,  as  in  this  case,  additional  movement  must  be  superimposed  on  the 
mesh  to  prevent  crossover  and/or  “piling  up”  of  mesh  points  adjacent  to  the  moving 
surface.  The  vector  x  was  then  separated  into  two  vectors 

x=x^+% 

such  that  x^  was  grid  motion  due  to  the  adaption  and  x^  was  grid  motion  required  by 
boundary  movement.  The  vector  x,,  remains  in  the  flow  solver  portion  of  the  split 

algorithm,  since  the  flow  must  respond  to  boundary  movement  in  the  inertial  frame. 
This  cell  side  velocity  was  determined  by  approximating  a  Laplacian  distribution  over 
the  domain.  The  vector  x^  was  determined  in  the  adaption  step  and  was  used  in  the 

solution  update  step  as  before.  In  order  to  ensure  that  grid  crossover  does  not  occur, 
movement  was  restricted  to  a  dual  cell  in  the  standard  adaptive  stencil. 

Boundary  conditions  for  both  fixed  and  moving  solid  walls  require  that  ghost  cell 
values  of  the  primitive  variables  be  set  such  that  there  is  no  mass  flux  through  a  solid 
wall.  For  inviscid  walls,  the  slip  condition  requires  that  the  ghost  cell  velocity  tangent  to 
the  wall  be  equal  to  the  tangential  velocity  of  the  first  interior  cell.  The  normal 
component  of  the  velocity  in  the  ghost  cell  must  be  defined  such  that  the  interpolated 
normal  component  of  the  fluid  velocity  at  the  wall  is  equal  to  the  normal  component  of 
the  wall  speed.  For  an  ri=constant  surface,  this  condition  yields  the  following  first  order 
relations  for  the  ghost  cell  values  of  velocity: 
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Here,  the  subscript  (s)  refers  to  the  wall,  (S+V2)  refers  to  the  ghost  cell  value,  and 
(S-I/2)  refers  to  the  value  in  the  first  interior  cell.  Similar  relations  are  derived  for 
^=constant  boundaries. 

For  a  viscous  wall,  the  no  slip  condition  must  be  satisfied.  The  ghost  cell  values 
of  velocity  must  be  set  such  that  the  fluid  and  wall  velocity  are  exactly  equal.  These 
conditions  yield  the  first  order  relationships; 

u  ,  =2x-m  , 

2  2 

and 

V  ,  =2y-v  , 

^*2  ■'  2 

for  viscous  walls. 

Solid  walls  were  assumed  to  be  adiabatic  and  the  normal  component  of  the 
momentum  equation  was  used  to  obtain  pressure.  The  remaining  conditions  for  inflow, 
outflow,  and  grid  block  boundaries  were  identical  to  those  implemented  by  Ingram. 

A  final  change  was  necessary  to  accommodate  moving  boundaries.  The 
parametric  space  adaption  was  continued  although  a  modification  was  required.  After  a 
solver  time  step  was  completed,  the  physical  boundary  has  now  been  allowed 
(potentially)  to  move  to  a  new  location.  If  the  boundary  point  moves  into  the  interior  of 
the  original  domain,  then  a  transformation  exists  to  remap  the  point  back  into  the 
physical  space  after  the  adaption  procedure.  However,  if  the  point  moves  outside  the 
domain,  no  transformation  exists  to  allow  the  correct  remapping.  A  solution  to  this 
problem  was  to  reset  the  parametric  space  back  to  its  original  state  (i.e.,  having  integer 
coordinates)  and  then  define  a  new  transformation  at  the  beginning  of  every  call  to  the 
adapter.  The  adapter  then  iterates  the  grid  to  the  appropriate  level  of  convergence  in 
the  new  parametric  space  and  correctly  maps  the  solution  back  to  the  physical  space, 
preserving  the  time  accuracy  of  the  prescribed  boundary  motion. 

The  weight  function  used  in  the  center  of  mass  algorithm  serves  to  control  the 
amount  and  location  of  point  clustering.  More  recent  work  by  Laflin  and  McRae  also 
stresses  the  improvements  obtained  by  aligning  the  grid  with  the  flow  features.  Laflin 
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derives  a  weight  function  based  on  a  solution  interpolation  error  approximation 
(SIERRA)  which  serves  to  both  refine  and  align  the  grid  to  the  solution. 

During  the  course  of  this  work,  a  variety  of  weight  function  formulations  were 
examined.  In  evaluating  the  weight  functions,  it  is  important  to  first  appreciate  the 
effects  of  the  rapid  movement  of  the  feature.  In  a  temporally  evolving  solution,  if  the 
grid  is  heavily  refined  over  only  a  cell  or  two,  then  the  flow  feature  may  potentially  move 
into  a  less  refined  region  during  the  next  time  step  thereby  mitigating  the  benefits  of  the 
adaption  process.  Ideally,  the  grid  must  be  sufficiently  refined  over  the  entire  region 
traversed  by  the  flow  feature  through  a  time  step.  Since  the  grid  adapter  does  not  have 
access  to  future  solution  data,  the  grid  cannot  be  adapted  to  provide  the  exact 
refinement  required  over  the  forthcoming  time  step.  Considering  that  the  weight 
function  is  in  essence  a  measure  of  the  desired  mesh  density,  an  arbitrary  weight 
function  should  contain  local  extrema  in  the  regions  requiring  the  most  grid  point 
movement.  Localized  smoothing  of  this  function  would  diffuse  the  extrema,  thereby 
producing  a  smoother  variation  in  cell  volumes  and  providing  an  increase  in  the  radius 
of  clustering  that  is  proportional  to  the  amount  of  smoothing.  With  sufficient  smoothing, 
the  resulting  grid  will  generally  be  adequately  refined  over  the  domain  traversed  by  a 
flow  feature  throughout  the  next  time  step. 

Examination  of  weight  functions  based  on  first  and  second  derivatives  of  the 
primitive  fluid  variables  revealed  that  both  provided  an  excellent  locator  of  regions 
requiring  adaption  at  a  minimum  of  computational  effort.  First  derivative  weight 
functions,  however,  do  possess  the  drawback  of  overadapting  the  interior  of  expansion 
waves.  Second  derivative,  or  curvature  based  weight  functions  help  to  alleviate  this 
problem. 

Experimentation  with  Laflin’s  SIERRA  weight  function  has  shown  to  provide 
results  almost  identical  to  those  obtained  by  more  cost  effective  curvature  methods  in 
one  dimensional  problems.  For  two  dimensional  problems,  the  SIERRA  weight  function 
tends  to  result  in  less  grid  shear  for  certain  problems.  For  the  nearly  1-D  solutions 
presented  herein,  the  increases  in  performance  do  not  seem  to  offset  the  extra 
computational  expense  however.  For  this  reason,  the  weight  function  used  for  this  work 
was  based  on  second  derivatives  evaluated  in  the  computational  space  and  seemed  to 
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provide  a  balance  between  cost  and  performance.  SIERRA  will  be  used  for  more 
complex  flows.  Separate  weight  functions  were  computed  for  each  coordinate  direction 
to  allow  greater  control  over  the  adaption  process. 

Square  Wave  Validation  of  Grid  Convergence  Requirement 

A  simple  numerical  experiment  was  used  to  demonstrate  the  effects  of  grid 
convergence.  An  analytically  defined  square  wave  was  allowed  to  propagate  through  a 
long  narrow  channel,  similar  to  the  movement  of  a  shock  wave  in  a  shock  tube,  and 
the  ability  of  the  adaption  algorithm  to  track  the  wave  was  observed.  The  square  wave 
had  a  peak  value  of  100,  a  minimum  value  of  1,  and  traveled  in  the  positive  x-direction 
with  velocity  1  on  an  initial  grid  of  201  □  21  points.  Figure  1  shows  the  movement  of 
this  wave  through  the  domain  and  the  adapted  grids  at  four  instances  in  time.  In  this 
case,  the  adaption  algorithm  makes  20  explicit  iterations  per  time  step.  The 
discontinuities  was  resolved  well  initially  and  the  remainder  of  the  domain  was  left  with  a 
generous  number  of  points.  As  the  discontinuities  continue  traveled  down  the  channel, 
increasingly  more  grid  points  were  drawn  from  upstream  to  resolve  the  advancing  wave 
fronts  and  then  remained  where  unneeded  in  the  downstream  region.  This  process  can 
lead  to  the  near  total  evacuation  of  points  from  the  region  upstream  of  the 
discontinuities.  If  a  Navier-Stokes  solver  were  coupled  with  a  grid  that  behaved  in  this 
manner,  the  increasing  cell  size  in  the  region  of  the  discontinuities  could  lead  to  an 
unstable  solution.  Note  that  this  phenomenum  is  not  a  problem  (and  in  fact  may  be 
desirable)  for  steady  or  slowly  varying  solutions. 

The  observed  grid  results  from  the  fact  that  the  grid  equations  were  not 
converged  to  a  solution  that  was  truly  elliptic  in  nature.  The  solution  to  this  problem  was 
to  bring  the  grid  closer  to  convergence  each  time  the  adaptor  was  applied. 
Unfortunately,  increasing  grid  convergence  increased  the  cost  of  obtaining  the  solution. 

Figure  2  presents  the  same  square  wave,  but  this  time  the  adaption  algorithm 
used  200  implicit  iterations  per  time  step.  In  this  case,  the  grid  remains  finely  adapted 
in  the  discontinuous  regions  and  evenly  distributed  throughout  the  remainder  of  the 
domain.  By  allowing  the  grid  to  advance  to  a  sufficiently  converged  solution,  any  local 
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changes  in  the  grid  were  felt  by  distant  grid  points  and  the  entire  grid  was  adjusted  to 
accommodate  the  local  refinement. 

Diaphragm  Rupture  Without  Adaption 

After  installation  of  moving  boundary  capability  and  other  code  improvements, 
the  first  step  leading  to  coupling  with  the  solid  mechanics  code  was  to  couple  with 
controlled  diaphragm  motion.  The  diaphragm  shape  and  motion  are  not  realistic  but 
were  chosen  to  allow  assessment  of  expected  problems  with  grid  behavior  upon 
rupture. 

A  four  block  grid,  each  block  having  75x50  grid  points,  was  used  to  model  the 
domain.  The  shock  tube  had  length  10,  a  radius  of  .5,  and  the  diaphragm  was 
assumed  to  have  ruptured  initially  without  motion  and  to  have  an  initial  opening  radius 
of  .05.  The  diaphragm  was  .05  units  thick  at  the  tunnel  wall  and  tapered  down  to  a 
sharp  point  near  the  centerline.  The  diaphragm  was  initially  symmetric  and  centered 
about  the  line  x=0.  Two  grid  blocks  abut  the  diaphragm  and  the  remaining  two  blocks 
span  the  distance  between  the  centerline  and  the  tip  of  the  diaphragm.  The 
deformation  of  the  diaphragm  was  governed  by  a  linear  decrease  in  y*  from  the  tip 

velocity  to  zero  at  the  wall  and  a  cubic  decrease  in  from  the  tip  velocity  to  zero  at  the 

wall.  Initial  conditions  for  this  case  were  similar  to  those  of  a  standard  shock  tube 
problem.  The  pressure  ratio  across  the  diaphragm  was  10:1 ,  the  density  ratio  was  8:1 , 
and  the  flow  was  initially  at  rest. 

Figure  3  shows  the  time  history  of  the  opening  process  from  t=0  to  t=0.99  for  the 
case  of  a  diaphragm  tip  velocity  of  0.1  and  no  grid  adaption.  The  typical  features 
associated  with  a  1-D  shock  tube  were  observed,  but  here  the  waves  were  curved 
because  the  streamlines  converge  while  approaching  the  diaphragm  and  diverge  after 
passing  through  the  diaphragm.  Following  a  streamline,  fluid  from  the  driver  section  of 
the  tube  accelerates  to  supersonic  velocities  through  the  small  opening  in  the 
diaphragm  before  being  slowed  to  subsonic  velocities  by  a  shock  further  downstream. 
The  contact  discontinuity  formed  at  flow  initiation  then  occurs.  The  jet  of  fluid  is 
separated  from  the  quiescent  flow  into  which  it  is  entering  by  a  bow  shock.  As  the  jet 
enters  the  driven  section,  it  entrains  fluid  from  behind  the  diaphragm  resulting  in  the 


12 


formation  of  a  vortex  behind  the  tip.  This  vortex  is  then  convected  downstream  into  the 
driven  section.  The  last  frame  in  Figure  3  shows  that  the  jet  has  expanded  to  reach  the 
outer  wall  of  the  tube  and  the  bow  shock  is  now  nearly  planar.  The  flow  along  the  wall 
of  the  tube  develops  a  boundary  layer  which  the  unadapted  grid  does  not  resolve. 

Figure  4  shows  the  time  history  of  the  opening  process  from  t=0  to  t=0.33  for  the 
case  of  a  diaphragm  tip  velocity  of  0.3  and  no  grid  adaption.  At  the  end  of  this  run,  the 
diaphragm  deformation  was  equal  to  that  at  the  end  of  the  run  in  Figure  3.  The  same 
features  observed  for  the  case  with  tip  velocity  of  0.1  develop  for  this  case  with  the 
exception  of  the  vortex  roll-up  behind  the  tip.  Because  the  induced  fluid  velocity  behind 
the  diaphragm  nearly  matches  the  jet  velocity,  the  formation  of  the  tip  vortex  is 
inhibited.  Additionally,  the  run  was  stopped  to  ensure  a  diaphragm  displacement 
equivalent  to  that  of  the  previous  case,  therefore  the  jet  has  not  yet  expanded  to  the 
diameter  of  the  tube. 

Diaphragm  Rupture  With  Adaption 

The  adaptive  algorithm,  as  originally  formulated  for  fixed  boundary  problems, 
does  not  yet  provide  useful  adaption  in  the  vicinity  of  the  simulated  diaphragm  tip.  The 
grid  block  topology  in  this  region  leads  to  a  concave  grid  boundary,  even  in  parametric 
space.  This  leads  to  severe  problems  with  grid  line  crossover  near  the  time  an 
exponentially  decaying  grid  movement  restriction  function  centered  at  the  tip  results  in 
the  interim  solution  of  Figure  5.  Analysis  of  these  results  indicates  that  the  blocks  must 
be  modified  to  allow  discontinuous  grid  lines  across  the  boundary.  This  will  allow 
independent  adaption  in  each  block  and  alleviate  much  of  the  skewness  evident  in 
Figure  5. 

It  is  apparent  from  these  cases  that  changing  the  rate  at  which  the  diaphragm 
opens  affects  the  flow  structure.  It  will  therefore  be  critical  in  future  work  to  accurately 
determine  how  the  diaphragm  deforms  in  response  to  the  dynamic  loads  imposed  by 
the  fluid.  In  preparation  for  coupling  with  the  solid  mechanics  code,  Figure  6  shows  an 
initial  mesh  as  it  deforms  with  the  initial  deformation  of  the  diaphragm  geometry  in  use 
for  the  solids  code  development.  Also,  numerical  experiments  were  conducted  to 
examine  the  flow  structure  after  the  shock  wave  reflects  from  the  tunnel  wall.  An 


13 


example  of  the  complex  structure  of  the  resulting  shock-wave  boundary  layer  interaction 
as  resolved  by  the  adaptive  mesh  code  is  shown  in  Figure  7.  Note  the  jet  of  flow  near 
the  wall  entering  the  stationary  fluid  region  (  similar  to  that  shown  in  Weber  and  Oran). 
The  status  of  the  solid  mechanics  portion  of  the  work  follows. 

SOLID  MECHANICS  SIMULATION 

During  the  solid  mechanics  research,  plane  strain  and  axisymmetric  finite- 
element  programs  have  been  developed  for  the  investigation  of  dynamic  inelastic 
deformation  and  rupture  of  a  shock-tube  diaphragm.  Modifications  have  also  been 
applied  for  the  characterization  of  damage  progression  and  failure  evolution  pertaining 
to  diaphragm  rupture.  The  dynamic  constitutive  formulation  that  were  developed 
account  for  high  pressures,  strain-rate,  and  thermal  effects.  Specialized  numerical 
techniques  have  been  developed  that  allow  tracking  of  propagating  stress  waves  in 
structural  components.  This  is  significant,  because  present  results  indicate  that  crack 
nucleation  occurs  when  compressive  waves  at  high  pressures  reflect  off  the  rear  of  the 
diaphragm  as  tensile  waves  that  are  a  precursor  to  crack  initiation.  This  two  dimensional 
wave  reflection  was  used  as  a  ductile  failure  criterion  to  signify  crack  initiation  for  high 
pressure  loading  conditions.  Using  this  failure  criterion,  stresses  that  correspond  to  the 
fracture  strength  of  the  material  have  been  predicted.  The  generalized  inelastic  rate- 
dependent  phenomenological  plasticity  formulation  that  accounts  for  noncoaxial  effects 
such  as  thermal  softening  and  kinematic  and  isotropic  hardening  has  been  used,  with 
the  newly  developed  failure  criteria,  to  investigate  damage  progression  and  failure 
evolution  in  rupturing  diaphragms.  Unnotched,  notched,  and  cracked  geometries  have 
been  used  to  analyze  the  inelastic  strain  fields  that  could  be  associated  with  either  the 
gradual  opening  (petal  opening),  or  the  instantaneous  rupturing  of  the  diaphragm.  New 
methods  have  developed  to  track  the  failure  due  to  the  reflected  waves.  Results  from 
this  study  indicate  that  stress  fields,  stress  gradients,  and  plastic  strains  at  the  free-edge 
or  at  the  crack  or  notch  are  physically  large  enough,  at  the  higher  pressures,  to  result  in 
crack  nucleation,  diaphragm  rupture,  and  material  separation.  This  buildup  in  the 
opening  mode  stress  is  directly  related  to  the  reflected  tensile  waves. 
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Computational  Method 

The  current  deformation  state  is  obtained  by  updating  the  total  deformation 
tensor,  the  plastic  deformation-rate  tensor,  and  the  equation  of  state  as  a  function  of  the 
current  pressure  and  particle  velocities.  The  total  deformation  rate  tensor  and  the  total 
spin  tensor  are  obtained  by  an  explicit  finite-element  method.  The  equations  of  motion 
are  integrated,  by  the  central  difference  method,  to  obtain  the  nodal  accelerations,  the 
nodal  velocities  and  the  nodal  displacements.  The  plastic  deformation-rate  tensor  is 
obtained  by  the  derivation  of  an  initial  value  system  of  equations.  This  is  based  on  a 
method  developed  by  Zikry  (1994).  This  initial  value  system  of  equations  is  integrated 
by  a  combination  of  explicit  multistep  and  implicit  A-stable  methods.  A  family  of  A-stable 
methods  is  used  where  numerical  stiffness  of  the  initial  value  system  may  arise  due  to 
plastic  wave  speeds  at  the  interface.  The  stiffness  ratio  is  used  to  automatically  switch 
solution  methods  between  the  explicit  and  the  A-stable  methods.  The  validity  of  the 
present  algorithm  has  been  tested  by  simulating  the  axisymmetric  deformation  of 
several  plasticity  problems  (no  noncoaxial  effects)  for  which  known  numerical  solutions 
exist. 

An  unnotched,  undeformed  diaphragm,  with  a  pre-existing  central  through-crack 
(Fig.  8A),  was  subjected  to  pressures  ranging  from  0  to  104  psi.  The  diaphragm  was 
modeled  as  a  cantilevered  plate  with  a  free  end,  to  simulate  the  through-crack.  A 
uniform  pressure  was  applied  transversely,  along  one  side.  Two  rates  of  pressure 
application  were  examined:  an  instantaneous  load  application  and  a  ramp  load.  This 
was  done  to  determine  the  effects  of  loading  rate.  For  the  slower  rate  of  loading,  the 
pressure  was  applied  linearly  from  an  undeformed  state  (Fig.  8B).  To  gain  a  more 
detailed  understanding  of  the  deformation  modes  associated  with  these  meshes, 
contours  of  effective  plastic  strains  are  given  at  pressures  of  3000  psi,  6000  psi,  and 
7000  psi  (Figs.  9-11).  As  shown  in  Fig.  10,  at  a  pressure  of  6000  psi  (30  ms),  the 
effective  plastic  strain  has  attained  a  maximum  value  of  more  than  43%  near  the  fixed 
edge.  As  these  contours  also  indicate,  large  plastic-strain  gradients  have  developed  in 
a  small  region  near  the  fixed  edge.  Maximum  effective  plastic  strains  range  from 
approximately  0%  to  70%,  when  the  diaphragm  penetrates  the  plane  of  the  fixed  edge 
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(at  approximately  35  ms).  Results  for  the  impact  load  were  qualitatively  similar,  but  the 
diaphragm  had  penetrated  the  plane  of  the  fixed  edge  at  approximately  15  ms(Fig.  11). 

The  deformed  mesh  at  face  pressures  of  3000  psi,  6000  psi,  and  7000  psi  are 
also  shown  in  Figs.  9-11.  A  highly  refined  mesh  was  used  near  the  fixed  end,  in  the 
region  of  high  plastic  strain  gradients,  and  in  the  crack-tip  region.  These  refined 
meshes  enabled  us  to  resolve  the  high  stress  gradients  associated  with  crack-tip 
plasticity  and  high  pressures. 

Impact 

Based  on  the  newly  developed  solid  mechanics  code  for  structural  integrity,  the 
onset  of  failure  was  predicted  and  used  to  model  petal  opening.  The  solid  mechanics 
code  should  be  coupled  interactively  to  the  fluid  mechanics  code  for  fully  three- 
dimensional  simulations  and  modeling  of  material  response. 

FUTURE  WORK 

Fluid  Dynamic  Code 

Further  development  of  the  adaptive  grid  algorithm  will  be  needed  to  resolve  the 
complex  geometry  and  flow  movement  associated  with  diaphragm  opening.  This  will 
take  two  primary  directions: 

1 )  Development  of  non-continuous  grid  block  boundaries  so  that  adaption  can  be 
independent  between  blocks.  This  should  alleviate  the  severe  restrictions  on 
adaption  at  concave  physical  comers. 

2)  Development  of  a  hybrid  structured/unstructured  grid  approach  to  improve  ability  to 
couple  with  geometry  during  the  final  stages  of  the  opening  and  also  the  3-D  rupture 
geomeries. 

3)  Development  of  coupling  protocols  with  the  solid  mechanics  code. 

Solid  Mechanics  Code 

1)  Development  of  coupling  protocols  with  the  fluid  mechanics  code. 

2)  Extension  of  the  code  to  3-D.. 
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Figure  5 :  Adapted  Mesh  Diaphragm 
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Figure  6 :  Initial  Coupling  of  Moving 
Mesh  and  Solid  Mechanics  Deformation 


Note:  every  other  grid  line  plotted  In  y-directlon 


Figure  8A,B:  Loading  Distribution  and  Onset 
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Figure  9:  Contour  of  Effective  Strain  at  3000psi 
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Figure  10:  Contour  of  Effective  Strain  at  6000psi 


Figure  11:  Contour  of  Effective  Strain  at  7000psi 


